function T = slfld(X, g, varargin)
%SLFLD Trains the Fisher Linear Discriminant Transform
%
% [ Syntax ]
%   - T = slfld(X, nums, ...)
%   - T = slfld(X, ginds, ...)
%
% [ Arguments ]
%   - X:        The sample matrix
%   - nums:     The numbers of samples in classes
%   - ginds:    The cell array of sample indices for classes
%   - T:        The learned transform matrix
%
% [ Description ]
%   - T = slfld(X, nums, ...) learns the linear discriminant transform
%     by maximizing Fisher's criteria. 
%
%     Suppose there are n samples in a d-dimensional space from nc classes.
%     Then X should be a d x n matrix, with each column representing a
%     sample. The samples from the same class should be arranged together.
%     nums is a vector of length nc, nums(k) is the number of the samples
%     from the k-th class.
%
%     Given a sample x and the learned transform T, the discriminant
%     feature can be obtained by 
%         $ y = T' * (x - vmean) $
%     where vmean is the mean vector of the sample space. Please note that
%     T should be transposed before being multiplied.
%
%     For a set of samples in a sample matrix X, the corresponding feature
%     set can be obtained by
%         $ Y = T' * bsxfun(@minus, X, vmean) $
%
%     The implementation follows the standard two-stage procedure:
%       # get the whitening transform W by solving eigen-problem of Sw;
%       # learn the subspace projection P from W' * Sb * W;
%       # combine the transform as T = W * P.
%
%     You can further specify the following options to control the
%     learning.
%     \{:
%        - prepca:          Whether to perform PCA preceding the learning
%                           (The preceding PCA preserves all dimensions up
%                            to the rank of the sample set)
%                           default = false.
%
%        - weights:         The weights of the samples. 
%                           By default, it is [], which means all samples
%                           share the same weight 1. You can specify the
%                           weights as a vector of length n.
%
%        - whiten:          A cell array in form of {method} or {method, r}
%                           to specify how to do the whitening of Sw.
%                              - {'energy', r} solves W such that the ratio
%                                r of the original space spanned by Sw is
%                                preserved in the whitened space.
%                                {'energy'} => {'energy', 0.99}
%
%                              - {'truncate', r} solves W with the trailing
%                                eigenvectors corresponding to the
%                                eigenvalues less than r * max(eigvals)
%                                {'truncate'} => {'truncate', 1e-6}
%                               
%                              - {'reg', r}, solves W on the regularized Sw
%                                which equals Sw + r * max(eigvals) * I,
%                                instead of the original Sw, and preserves
%                                the whole space spanned by the regularized
%                                Sw. {'reg'} => {'reg', 1e-4}
%
%                              - {'bound', r} solves W on a regularized Sw,
%                                which is formed by forcing all trailing 
%                                eigenvalues that are less than the lower
%                                bound lb = r * max(eigvals) to lb.
%                                {'bound'} => {'bound', 1e-4}
%
%                              - {'avgtr', r} solves W on a regularized Sw,
%                                which is formed by setting all trailing
%                                eigenvalues that are less than the
%                                threshold r * max(eigvals) to their
%                                average. {'avgtr'} => {'avgtr', 1e-4}
%
%                           By default, it is set to {'reg', 1e-4}.
%
%        - preserve:        A cell array in form of {cri} or {cri, v} to specify 
%                           how many features should be preserved in the
%                           output transform.
%
%                           Here are the available criteria.
%                           \{:
%                             - rank:   preserve all eigenvalues not
%                                       smaller than eps(max(eigvals))
%                             - num:    preserve all the v leading
%                                       eigenvalues
%                             - ratio:  preserve all eigenvalues not
%                                       smaller than (v * max(eigvals)
%                             - energy: preserve energy of ratio not less
%                                       than v.
%                           \:}
%                           By default, the option is set to {'rank'}.
%     \:}
%
%   - T = slfld(X, ginds, ...) uses the indices of different classes to
%     specify which samples are from which classes. 
%
%     ginds is a cell array with nc elements. ginds{k} is the array of the
%     indices of the samples from the k-th class.
%
% [ History ]
%   - Created by Dahua Lin, on Aug 16, 2006
%   - Modified by Dahua Lin, on Jul 23, 2006
%       - rewrite the whole program
%       - clarify the help
%       - remove the support of pre-computed Sb and Sw.
%       - new whiten and preserve options.
%

%% parse and verify input arguments

error(nargchk(2, inf, nargin));

% X
assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slfld:invalidarg', ...
    'X should be a 2D numeric matrix.');
n = size(X, 2);

% g (nums | ginds)
assert((isnumeric(g) || iscell(g)) && isvector(g), 'sltoolbox:slfld:invalidarg', ...
    'The 2nd argument can be a numeric vector as nums or a cell vector as ginds.');
nc = length(g);
if ~isequal(size(g), [nc 1])
    g = reshape(g, [nc 1]);
end
if isnumeric(g)
    assert(sum(g) == n, 'sltoolbox:slfld:nummismatch', ...
        'The number of samples given in nums is not equal to that given in X.');    
    eps = cumsum(g);
    sps = [1; eps(1:end-1)+1];
    ginds = arrayfun(@(sp, ep) sp:ep, sps, eps, 'UniformOutput', false);
else
    ginds = g;
end

% options

opts = struct( ...
    'prepca', false, ...
    'weights', [], ...
    'whiten', {{'reg', 1e-4}}, ...
    'preserve', {{'rank'}});
if ~isempty(varargin)
    opts = setopts(opts, varargin{:});
end

prepca = opts.prepca;
assert(islogical(prepca) && isscalar(prepca), 'sltoolbox:slfld:invalidopt', ...
    'The option prepca should be a logical scalar.');

weights = opts.weights;
if ~isempty(weights)
    assert(isnumeric(weights) && isvector(weights) && length(weights) == n, ...
        'sltoolbox:slfld:invalidopt', ...
        'The option weights should be a numeric vector of length n.');
end

whiten = opts.whiten;
if ~isempty(whiten)
    assert(iscell(whiten) && (numel(whiten) == 1 || numel(whiten) == 2), ...
        'sltoolbox:slfld:invalidopt', ...
        'The option whiten should be a cell vector of length 1 or 2.');
end
f_whiten_evs = fld_whiten_fun(whiten{:});

preserve = opts.preserve;
if ~isempty(preserve)
    assert(iscell(preserve) && (numel(preserve) == 1 || numel(preserve) == 2), ...
        'sltoolbox:slfld:invalidopt', ...
        'The option preserve should be a cell vector of length 1 or 2');
end
f_preserve = geteigpreserve(preserve{:});


%% main

% preparation

% prepca (if specified)
if prepca
    SPCA = slpca(X, 'weights', weights);
    X = slpcafea(SPCA, X);
    clear SPCA;
end

% compute the class-specific means and weights
if isempty(weights)
    cmeans = cellfun(@(u) slmean(X(:, u)), ginds, 'UniformOutput', false);
    cweights = cellfun(@numel, ginds);
else
    cmeans = cellfun(@(u) slmean(X(:, u), weights(u)), 'UniformOutput', false);
    cweights = cellfun(@(u) sum(weights(u)), ginds);
end
cmeans = [cmeans{:}];

% step 1: compute the Sw-whitening transform W

W = fld_whiten(X, weights, ginds, cmeans, f_whiten_evs);

% step 2: whiten Sb and solve optimal subspace projection

whiten_cmeans = W' * cmeans;
wSb = slcov(whiten_cmeans, cweights);
[evals, P] = slsymeig(wSb); 

% preserve 
if ~isempty(f_preserve)
    inds = f_preserve(evals);
    P = P(:, inds);
end

% step 3: combine the two-stage transformation

T = W * P;



%% The functions for whitening

function W = fld_whiten(X, w, ginds, cmeans, f_whiten_evs)

mrk = min(size(X,1), size(X,2)-1);
[D, wd] = make_withinclass_diffvecs(X, ginds, cmeans, w);
Sw = slcov(D, wd, 0);
clear D wd;

[evals, V] = slsymeig(Sw);
if size(V,2) > mrk
    evals = evals(mrk);
    V = V(:, mrk);
end
evals(evals < 0) = 0;

if ~isempty(f_whiten_evs)
    evals = f_whiten_evs(evals);
    if numel(evals) < size(V, 2)
        V = V(:, 1:numel(evals));
    end
end

evals = evals(:)';
W = bsxfun(@times, V, 1 ./ evals);



function f = fld_whiten_fun(method, r)

assert(ischar(method), 'sltoolbox:slfld:invalidopt', ...
    'The first argument for option whiten should be a string of the method name.');

switch method
    case 'energy'
        if nargin < 2
            r = 0.99;
        else
            assert(isnumeric(r) && isscalar(r) && isreal(r) && r > 0 && r < 1, ...
                'sltoobox:slfld:invalidopt', ...
                'The parameter for energy should be a real scalar ranging in (0, 1).');
        end
        f = @(x) x(1:find(cumsum(x) >= sum(x) * r, 1));
        
    case 'truncate'
        if nargin < 2
            r = 1e-6;
        else
            assert(isnumeric(r) && isscalar(r) && isreal(r) && r > 0 && r < 1, ...
                'sltoobox:slfld:invalidopt', ...
                'The parameter for truncate should be a real scalar ranging in (0, 1).');
        end
        f = @(x) x(1:find(x >= x(1) * r, 1, 'last'));
            
    case 'reg'
        if nargin < 2
            r = 1e-4;
        else
            assert(isnumeric(r) && isscalar(r) && r > 0, ...
                'sltoobox:slfld:invalidopt', ...
                'The parameter for reg should be a positive real scalar.');
        end
        f = @(x) x + r * x(1);
        
    case 'bound'
        if nargin < 2
            r = 1e-4;
        else
            assert(isnumeric(r) && isscalar(r) && r > 0 && r < 1, ...
                'sltoobox:slfld:invalidopt', ...
                'The parameter for bound should be a real scalar ranging in (0, 1).');
        end
        f = @(x) max(x, r * x(1));
        
    case 'avgtr'
        if nargin < 2 
            r = 1e-4;
        else
            assert(isnumeric(r) && isscalar(r) && r > 0 && r < 1, ...
                'sltoobox:slfld:invalidopt', ...
                'The parameter for energy should be a real scalar ranging in (0, 1).');
        end
        f = @(x) fld_whiten_avgtr(x, r);
        
    otherwise
        error('sltoolbox:slfld:unknownmethod', ...
            'The whitening method %s is unknown', method);
end


function xr = fld_whiten_avgtr(x, r)

p = find(x < r * x(1));
xr = x;
if ~isempty(p)
    xr(p) = sum(xr(p)) / numel(p);
end

